Authors: MJ Salazar-Nicholls1,3, KD Escobar2 & KM Warkentin3 1. Pontificia Universidad Católica, Ecuador; 2. Western Connecticut State University; 3. Boston Affiliations: University – majonicholls2909@gmail.com – Warkentin Lab
Agalychnis callidryas
wasp attacking Agalychnis callidryas clutch
Hatching early allows embryos to escape threats to eggs, but increases risks to larvae. Red-eyed treefrogs, Agalychnis callidryas, hatch by rapidly releasing enzymes to digest a small hole in their membrane, then squeezing out aided by turgor pressure. Displacement from the initial hole can occur spontaneously and in predator attacks, complicating hatching by capsular collapse as fluid escapes.
Assess developmental changes in ability to recover from hatching complications
| Flooding | Jiggling |
|---|---|
| we submerged younger eggs in hypoxic water (N = 25, 27 at ages 3 and 4 d) | we used a blunt probe to move older eggs on clutches (N = 64, 51 at ages 4 and 5 d) |
R1: rupture one
R2: rupture two
This is an analysis of the data used for a poster presentation at SICB 2017, which now is part of a manuscript. The statistical analyses presented in the poster were made using JMP-data analysis sorfware
#install.packages("MASS")
#install.packages("coin")
library(curl)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library(EnvStats)
##
## Attaching package: 'EnvStats'
## The following objects are masked from 'package:stats':
##
## predict, predict.lm
## The following object is masked from 'package:base':
##
## print.default
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:EnvStats':
##
## boxcox
## The following object is masked from 'package:plotly':
##
## select
## The following object is masked from 'package:dplyr':
##
## select
library(coin)
## Loading required package: survival
f <- curl("https://raw.githubusercontent.com/mjsalnic/mjsalnic-data-replication-assignment/master/HC2016R.csv")
HC2016R <- read.csv(f, header = TRUE, sep = ";",dec=",", stringsAsFactors = FALSE)
head(HC2016R)
## Treatment ClutchN Age Stage comst groupst DecimaClocklHours
## 1 Flooding 462 3 2 2-3 2-3 20.23
## 2 Flooding 470 3 3 2-3 2-3 21.02
## 3 Flooding 386 4 4 4-5 4-7 0.15
## 4 Flooding 397 4 5 4-5 4-7 0.82
## 5 Flooding 102 4 5 4-5 4-7 10.85
## 6 Flooding 187 4 4 4-5 4-7 16.12
## AccumshakingtimeRup1Fra AccumShakingR1sec NumPosChanADBefore2ndHole
## 1 4 0.13 3
## 2 290 9.67 7
## 3 30 1.00 1
## 4 395 13.17 1
## 5 105 3.50 17
## 6 260 8.67 1
## AccumADshakingR AccumshakingtimeRup2Fra AccumShakingR2sec DisEndRup2Fr
## 1 NA NA NA NA
## 2 NA 0 0 NA
## 3 NA 0 0 NA
## 4 NA NA NA NA
## 5 NA 0 0 NA
## 6 NA NA NA NA
## DisEndHatchingStartFr CompDurFr CompressionSec CompressionSeen
## 1 1407 NA NA 0
## 2 4838 NA NA 0
## 3 916 6 0.20 1
## 4 5182 NA NA 0
## 5 6511 7 0.23 1
## 6 2222 NA NA 0
## EXITDURATIONFr EXIT.DURATION.StarttoFullyOutsec NumHoles
## 1 2442 81.4000000 1
## 2 1318 43.9333333 1
## 3 342 11.4000000 1
## 4 4 0.1333333 1
## 5 3995 133.1666667 1
## 6 157 5.2333333 1
## FirstSecondHoleExit Rescue X X.1
## 1 First 0
## 2 First 0
## 3 First 0
## 4 First 0
## 5 First 0
## 6 First 0
attach(HC2016R)
Transforming numerical variables to factorial ### Exit duration
HC2016R$Stage<- factor(HC2016R$Stage)
HC2016R$Age<- factor(HC2016R$Age)
HC2016R$NumHoles<- factor(HC2016R$NumHoles)
HC2016R$Rescue<- factor(HC2016R$Rescue)
str(HC2016R)
## 'data.frame': 167 obs. of 25 variables:
## $ Treatment : chr "Flooding" "Flooding" "Flooding" "Flooding" ...
## $ ClutchN : int 462 470 386 397 102 187 272 268 239 244 ...
## $ Age : Factor w/ 3 levels "3","4","5": 1 1 2 2 2 2 2 2 2 2 ...
## $ Stage : Factor w/ 6 levels "2","3","4","5",..: 1 2 3 4 4 3 4 3 5 5 ...
## $ comst : chr "2-3" "2-3" "4-5" "4-5" ...
## $ groupst : chr "2-3" "2-3" "4-7" "4-7" ...
## $ DecimaClocklHours : num 20.23 21.02 0.15 0.82 10.85 ...
## $ AccumshakingtimeRup1Fra : int 4 290 30 395 105 260 285 0 307 170 ...
## $ AccumShakingR1sec : num 0.13 9.67 1 13.17 3.5 ...
## $ NumPosChanADBefore2ndHole : int 3 7 1 1 17 1 1 14 16 1 ...
## $ AccumADshakingR : int NA NA NA NA NA NA NA NA NA NA ...
## $ AccumshakingtimeRup2Fra : int NA 0 0 NA 0 NA NA NA 385 NA ...
## $ AccumShakingR2sec : num NA 0 0 NA 0 ...
## $ DisEndRup2Fr : int NA NA NA NA NA NA NA NA NA NA ...
## $ DisEndHatchingStartFr : int 1407 4838 916 5182 6511 2222 213 8227 5455 187 ...
## $ CompDurFr : int NA NA 6 NA 7 NA NA NA NA NA ...
## $ CompressionSec : num NA NA 0.2 NA 0.23 NA NA NA NA NA ...
## $ CompressionSeen : int 0 0 1 0 1 0 0 0 0 0 ...
## $ EXITDURATIONFr : int 2442 1318 342 4 3995 157 96 NA NA 147 ...
## $ EXIT.DURATION.StarttoFullyOutsec: num 81.4 43.933 11.4 0.133 133.167 ...
## $ NumHoles : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 2 1 ...
## $ FirstSecondHoleExit : chr "First" "First" "First" "First" ...
## $ Rescue : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ X : chr "" "" "" "" ...
## $ X.1 : chr "" "" "" "" ...
let’s subset the data into 2 outcome groups: success_hat (hatched embryos after displacement) and failed_hat(embryos unable to hatch after displacement)
success_hat<-subset(HC2016R, Rescue==0, na.rm=T)
failed_hat<-subset(HC2016R, Rescue==1, na.rm=T)
let’s subset the data into 2 treatment groups: flooding and jiggling
flooding<-subset (HC2016R, Treatment=="Flooding", na.rm=T)
jiggling<-subset(HC2016R, Treatment=="Jiggling", na.rm=T)
Explore if there are differences in the outcome between treatments ## Results
There are not differences between treatments
totals=table(HC2016R$Stage,HC2016R$Rescue)
mosaicplot(totals, main = "Escape-hatching success after complications", xlab = "Embryonic developmental stage", ylab = "proportion of embryos", col=c("light green", "black"))
#where 0 means no rescue needed, embryo able to hatch by itself and 1 means rescue after displacement experiment
1= Failed to hatch 0= Succesfully hatched embryos
Stage: P < 0.0001
tab1<-table(HC2016R$Stage, HC2016R$Rescue)
prop.table(tab1)*100
##
## 0 1
## 2 2.395210 4.191617
## 3 5.988024 1.796407
## 4 11.377246 0.000000
## 5 20.359281 1.796407
## 6 20.958084 0.000000
## 7 31.137725 0.000000
chisq.test(tab1)
## Warning in chisq.test(tab1): Chi-squared approximation may be incorrect
##
## Pearson's Chi-squared test
##
## data: tab1
## X-squared = 60.99, df = 5, p-value = 7.588e-12
fisher.test(tab1)
##
## Fisher's Exact Test for Count Data
##
## data: tab1
## p-value = 9.464e-09
## alternative hypothesis: two.sided
tab<-table(HC2016R$Stage, HC2016R$NumHoles)
prop.table(tab)*100
##
## 1 2
## 2 4.790419 1.796407
## 3 1.796407 5.988024
## 4 2.395210 8.982036
## 5 7.185629 14.970060
## 6 8.982036 11.976048
## 7 7.784431 23.353293
chisq.test(tab)
## Warning in chisq.test(tab): Chi-squared approximation may be incorrect
##
## Pearson's Chi-squared test
##
## data: tab
## X-squared = 12.719, df = 5, p-value = 0.02616
fisher.test(tab)
##
## Fisher's Exact Test for Count Data
##
## data: tab
## p-value = 0.03435
## alternative hypothesis: two.sided
attach(flooding)
## The following objects are masked from HC2016R:
##
## AccumADshakingR, AccumShakingR1sec, AccumShakingR2sec,
## AccumshakingtimeRup1Fra, AccumshakingtimeRup2Fra, Age,
## ClutchN, CompDurFr, CompressionSec, CompressionSeen, comst,
## DecimaClocklHours, DisEndHatchingStartFr, DisEndRup2Fr,
## EXIT.DURATION.StarttoFullyOutsec, EXITDURATIONFr,
## FirstSecondHoleExit, groupst, NumHoles,
## NumPosChanADBefore2ndHole, Rescue, Stage, Treatment, X, X.1
totals=table(flooding$Stage,flooding$NumHoles)
mosaicplot(totals, main = "Multiple membrane ruptures", xlab = "Embryonic developmental stage", ylab = "Proportion of embryos", col=c("light green", "black"))
p-value = 0.009272
2.1 Contigency table- Multiple membrane ruptures
tab2<-table(flooding$Stage, flooding$NumHoles)
fisher.test(tab2)
##
## Fisher's Exact Test for Count Data
##
## data: tab2
## p-value = 0.009272
## alternative hypothesis: two.sided
totals=table(jiggling$Stage,jiggling$NumHoles)
mosaicplot(totals, main = "Multiple membrane ruptures", xlab = "Embryonic developmental stage", ylab = "Proportion of embryos", col=c("light green", "black"))
p-value = 0.1508
3.1 Contigency table- Multiple membrane ruptures
tab3<-table(jiggling$Stage, jiggling$NumHoles)
fisher.test(tab3)
##
## Fisher's Exact Test for Count Data
##
## data: tab3
## p-value = 0.1508
## alternative hypothesis: two.sided
summary(tab3)
## Number of cases in table: 115
## Number of factors: 2
## Test for independence of all factors:
## Chisq = NaN, df = 5, p-value = NA
## Chi-squared approximation may be incorrect
Should we use parametric or non-parametric test?
par(mfrow=c(3,2))
hist(HC2016R$NumPosChanADBefore2ndHole[HC2016R$Stage=="2"])
hist(HC2016R$NumPosChanADBefore2ndHole[HC2016R$Stage=="3"])
hist(HC2016R$NumPosChanADBefore2ndHole[HC2016R$Stage=="4"])
hist(HC2016R$NumPosChanADBefore2ndHole[HC2016R$Stage=="5"])
hist(HC2016R$NumPosChanADBefore2ndHole[HC2016R$Stage=="6"])
hist(HC2016R$NumPosChanADBefore2ndHole[HC2016R$Stage=="7"])
par(mfrow=c(1,2))
boxplot(NumPosChanADBefore2ndHole~Stage,success_hat)
boxplot(NumPosChanADBefore2ndHole~Stage,failed_hat)
Subsetting data frame: to extract only values from the earliest stages where embryos got more trouble making a second hole to hatch
young_outcome<- subset(HC2016R, comst="2-3", select = NumPosChanADBefore2ndHole:Rescue)
str(young_outcome)
## 'data.frame': 167 obs. of 14 variables:
## $ NumPosChanADBefore2ndHole : int 3 7 1 1 17 1 1 14 16 1 ...
## $ AccumADshakingR : int NA NA NA NA NA NA NA NA NA NA ...
## $ AccumshakingtimeRup2Fra : int NA 0 0 NA 0 NA NA NA 385 NA ...
## $ AccumShakingR2sec : num NA 0 0 NA 0 ...
## $ DisEndRup2Fr : int NA NA NA NA NA NA NA NA NA NA ...
## $ DisEndHatchingStartFr : int 1407 4838 916 5182 6511 2222 213 8227 5455 187 ...
## $ CompDurFr : int NA NA 6 NA 7 NA NA NA NA NA ...
## $ CompressionSec : num NA NA 0.2 NA 0.23 NA NA NA NA NA ...
## $ CompressionSeen : int 0 0 1 0 1 0 0 0 0 0 ...
## $ EXITDURATIONFr : int 2442 1318 342 4 3995 157 96 NA NA 147 ...
## $ EXIT.DURATION.StarttoFullyOutsec: num 81.4 43.933 11.4 0.133 133.167 ...
## $ NumHoles : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 2 1 ...
## $ FirstSecondHoleExit : chr "First" "First" "First" "First" ...
## $ Rescue : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
par(mfrow=c(1,2))
hist(young_outcome$NumPosChanADBefore2ndHole[young_outcome$Rescue=="1"])
hist(young_outcome$NumPosChanADBefore2ndHole[young_outcome$Rescue=="0"])
boxplot(data = young_outcome, NumPosChanADBefore2ndHole ~ Rescue, col = c("burlywood2",
"lightpink1"))
wilcox.test(NumPosChanADBefore2ndHole~Rescue, data = young_outcome, correct=FALSE)
##
## Wilcoxon rank sum test
##
## data: NumPosChanADBefore2ndHole by Rescue
## W = 160, p-value = 8.265e-08
## alternative hypothesis: true location shift is not equal to 0
p-value = < 0.0001 N= 10,14
boxplot(data = success_hat, NumPosChanADBefore2ndHole ~ comst, col = c("burlywood2",
"lightpink1", "purple"))
kruskal.test(NumPosChanADBefore2ndHole ~ comst,
data =success_hat)
##
## Kruskal-Wallis rank sum test
##
## data: NumPosChanADBefore2ndHole by comst
## Kruskal-Wallis chi-squared = 2.7825, df = 2, p-value = 0.2488
p-value = 0.2488 N = 14, 53, 84
Exit from egg: fist or second rupture
boxplot(data = success_hat, NumPosChanADBefore2ndHole ~ NumHoles, col = c("burlywood2",
"lightpink1"))
wilcox.test(NumPosChanADBefore2ndHole~NumHoles, data = success_hat, correct=FALSE)
##
## Wilcoxon rank sum test
##
## data: NumPosChanADBefore2ndHole by NumHoles
## W = 3450.5, p-value = 3.33e-06
## alternative hypothesis: true location shift is not equal to 0
P < 0.0001 N = 41, 110
Exit duration
#good one
attach(success_hat)
## The following objects are masked from flooding:
##
## AccumADshakingR, AccumShakingR1sec, AccumShakingR2sec,
## AccumshakingtimeRup1Fra, AccumshakingtimeRup2Fra, Age,
## ClutchN, CompDurFr, CompressionSec, CompressionSeen, comst,
## DecimaClocklHours, DisEndHatchingStartFr, DisEndRup2Fr,
## EXIT.DURATION.StarttoFullyOutsec, EXITDURATIONFr,
## FirstSecondHoleExit, groupst, NumHoles,
## NumPosChanADBefore2ndHole, Rescue, Stage, Treatment, X, X.1
## The following objects are masked from HC2016R:
##
## AccumADshakingR, AccumShakingR1sec, AccumShakingR2sec,
## AccumshakingtimeRup1Fra, AccumshakingtimeRup2Fra, Age,
## ClutchN, CompDurFr, CompressionSec, CompressionSeen, comst,
## DecimaClocklHours, DisEndHatchingStartFr, DisEndRup2Fr,
## EXIT.DURATION.StarttoFullyOutsec, EXITDURATIONFr,
## FirstSecondHoleExit, groupst, NumHoles,
## NumPosChanADBefore2ndHole, Rescue, Stage, Treatment, X, X.1
#good one
p <- ggplot(data=success_hat, aes(x=Stage, y=EXIT.DURATION.StarttoFullyOutsec, fill=Stage))+stat_boxplot(geom='errorbar', linetype=1, width=0.5)+ #whiskers
geom_boxplot(outlier.shape=1)
p <- p + geom_boxplot()+ stat_summary(fun.y = mean, geom = "point", shape=5, size=5, fill="white")
p <- p + theme(axis.text.x=element_text(angle=0))
p <- p + ylab("Exit duration (sec)")
p <- p + stat_n_text(y.pos = 170)
p <- p + theme_classic()+ scale_fill_manual(breaks = c("2", "3", "4", "5", "6", "7"),
values=c("lightskyblue", "lightskyblue", "mediumaquamarine","mediumturquoise" ,"mediumturquoise","olivedrab2"))
p
## Warning: Removed 42 rows containing non-finite values (stat_boxplot).
## Warning: Removed 42 rows containing non-finite values (stat_boxplot).
## Warning: Removed 42 rows containing non-finite values (stat_boxplot).
## Warning: Removed 42 rows containing non-finite values (stat_summary).
## Warning: Removed 42 rows containing non-finite values (stat_n_text).
Body Compression
totals1=table(success_hat$comst,success_hat$CompressionSeen)
mosaicplot(totals1, main = "Escape-hatching success after complications", xlab = "Embryonic developmental stage", ylab = "proportion of embryos", col=c("light green", "pink"))
p <- ggplot(data=success_hat, aes(x=groupst, y= CompressionSec, fill=groupst))+stat_boxplot(geom='errorbar', linetype=1, width=0.5)+ #whiskers
geom_boxplot(outlier.shape=1)
p <- p + geom_boxplot()+ stat_summary(fun.y = mean, geom = "point", shape=5, size=5, fill="white")
p <- p + theme(axis.text.x=element_text(angle=0))
p <- p + ylab("Compression Duration (sec)") + xlab("Developmental stage")
p <- p + stat_n_text(y.pos = 70)
p <- p + theme_classic()+ scale_fill_manual(breaks = c("2,3", "4,7"),
values=c("lightskyblue", "mediumaquamarine"))
p
## Warning: Removed 124 rows containing non-finite values (stat_boxplot).
## Warning: Removed 124 rows containing non-finite values (stat_boxplot).
## Warning: Removed 124 rows containing non-finite values (stat_boxplot).
## Warning: Removed 124 rows containing non-finite values (stat_summary).
## Warning: Removed 124 rows containing non-finite values (stat_n_text).
detach(flooding)
attach(success_hat)
## The following objects are masked from success_hat (pos = 3):
##
## AccumADshakingR, AccumShakingR1sec, AccumShakingR2sec,
## AccumshakingtimeRup1Fra, AccumshakingtimeRup2Fra, Age,
## ClutchN, CompDurFr, CompressionSec, CompressionSeen, comst,
## DecimaClocklHours, DisEndHatchingStartFr, DisEndRup2Fr,
## EXIT.DURATION.StarttoFullyOutsec, EXITDURATIONFr,
## FirstSecondHoleExit, groupst, NumHoles,
## NumPosChanADBefore2ndHole, Rescue, Stage, Treatment, X, X.1
## The following objects are masked from HC2016R:
##
## AccumADshakingR, AccumShakingR1sec, AccumShakingR2sec,
## AccumshakingtimeRup1Fra, AccumshakingtimeRup2Fra, Age,
## ClutchN, CompDurFr, CompressionSec, CompressionSeen, comst,
## DecimaClocklHours, DisEndHatchingStartFr, DisEndRup2Fr,
## EXIT.DURATION.StarttoFullyOutsec, EXITDURATIONFr,
## FirstSecondHoleExit, groupst, NumHoles,
## NumPosChanADBefore2ndHole, Rescue, Stage, Treatment, X, X.1
x <- c(Stage)
x<- factor(x)
levels(x)<- list("St2-3"=c("2","3"), "st 4-7"=c("4", "5", "6", "7"))
x
## [1] <NA> St2-3 St2-3 st 4-7 st 4-7 St2-3 st 4-7 St2-3 st 4-7 st 4-7
## [11] st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7
## [21] st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 St2-3
## [31] st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7
## [41] st 4-7 st 4-7 st 4-7 st 4-7 <NA> St2-3 St2-3 St2-3 <NA> <NA>
## [51] St2-3 St2-3 St2-3 St2-3 St2-3 St2-3 st 4-7 st 4-7 St2-3 St2-3
## [61] St2-3 St2-3 St2-3 St2-3 St2-3 St2-3 St2-3 St2-3 St2-3 St2-3
## [71] st 4-7 St2-3 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 St2-3
## [81] st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 St2-3 st 4-7 st 4-7 st 4-7
## [91] st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7
## [101] st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7
## [111] st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7
## [121] st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7
## [131] st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7
## [141] st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7 st 4-7
## [151] st 4-7 st 4-7 st 4-7 st 4-7
## Levels: St2-3 st 4-7
boxplot(CompressionSec ~ x, data=success_hat)
p-value = 0.02087
wilcox.test(CompressionSec ~groupst, data=success_hat)
## Warning in wilcox.test.default(x = c(53.27, 7.5, 19.4, 15.4, 9.53, 0.03, :
## cannot compute exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: CompressionSec by groupst
## W = 128, p-value = 0.02087
## alternative hypothesis: true location shift is not equal to 0
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)
## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.